Statistical inferences under step stress partially accelerated life testing based on multiple censoring approaches using simulated and real-life engineering data

Evaluating the lifespan distribution of highly reliable commodities under regular use is exceedingly difficult, time consuming, and extremely expensive. As a result of its ability to provide more failure data faster and at a lower experimental cost, accelerated life testing has become increasingly important in life testing studies. In this article, we concentrate on parametric inference for step stress partially life testing utilizing multiple censored data based on the Tampered Random Variable model. Under normal stress circumstances, the lifespan of the experimental units is assumed to follow the Nadarajah–Haghighi distribution, with and being the shape and scale parameters, respectively. Maximum likelihood estimates for model parameters and acceleration factor are developed using multiple censored data. We build asymptotic confidence intervals for the unknown parameters using the observed Fisher information matrix. To demonstrate the applicability of the different methodologies, an actual data set based on the timings of subsequent failures of consecutive air conditioning system failures for each member of a Boeing 720 jet aircraft fleet is investigated. Finally, thorough simulation studies utilizing various censoring strategies are performed to evaluate the estimate procedure performance. Several sample sizes were studied in order to investigate the finite sample features of the considered estimators. According to our numerical findings, the values of mean squared errors and average asymptotic confidence intervals lengths drop as sample size increases. Furthermore, when the censoring level is reduced, the considered estimates of the parameters approach their genuine values.


Scientific Reports
| (2023) 13:12452 | https://doi.org/10.1038/s41598-023-39170-x www.nature.com/scientificreports/ the parameters α, β and θ are then addressed. To demonstrate the applicability of the various methodologies, an actual data set based on air conditioning system failure times for each member of a Boeing 720 jet aircraft fleet is evaluated in "Real engineering application" section. A thorough numerical analysis is performed, illustrating the positive behavior of the derived estimates over a wide range of sample sizes in "Simulation study" section.

Modeling SSPALT with MCS
Let Y is a nonnegative random variable distributed according to NH distribution with scale parameter β and shape parameter α , denoted as NH ( α, β ), then its probability density function (PrDF), cumulative distribution function (CDF) and the survival function (SF) are as follows: Different shapes of PrDF, CDF, and SF that were created using different input values of parameters are displayed in Fig. 1.
As a special instance of the NH distribution, when α = 1 , an exponential distribution can be produced. It offers closed versions of survival and hazard rate functions, such as that of the Weibull distribution, which makes it an excellent alternative for lifetime data investigators.
In SSPALT, all testing items are first allocated to be tested under ordinary usage settings until a pre-set stress change time τ , following which any remaining survivals that have not failed by time τ are moved to be tested under accelerated conditions. The effect of stress transition from normal to accelerated condition may be explained by multiplying the remaining lifetime by the inverse of the acceleration factor. The following is a theoretical calculation of the item's total lifetime T under accelerated conditions: where Y denotes the item's lifespan under ordinary usage settings and θ > 1 denotes the AF, which is often depending on applied stress. We can now explain the PDF, CDF, and RF of total life T of the objects using the model provided in Eq. (4) as follows: (1) f y; α, β = αβ 1 + βy α−1 exp 1 − 1 + βy α , y > 0, α, β > 0 (2) F y; α, β = 1 − exp 1 − 1 + βy α , y > 0, α, β > 0 (3) R y = exp 1 − 1 + βy α , y > 0, α, β > 0.  www.nature.com/scientificreports/ Assume we're dealing with an SSPALT based on MCS. Assume the test is based on only two stress levels, X 1 representing regular working conditions and X 2 representing accelerated conditions, where X 2 is larger than X 1 . This life testing consists of n objects that are both identical and independent in nature. Under each of the stress levels X 1 and X 2 , at least one failure should occur. The failure times of the items at each of the stress levels X 1 and X 2 are determined by the NH ( α, β ) presented by (1). All items in the sample of size n are now allocated to be tested at stress level X 1 with a known number r 1 of failures and a corresponding number m 1 of multiply censored items till time τ . The test will now be continued by testing all of the items from n that have not failed or been censored up to time τ at the stress X 2 with the pre-requisite r 2 number of failures and the associated number of multiply censored items m 2 until all of the test items have failed or been censored.

Inferences under SSPALT with MCS
In this section, the MLE method is utilized to estimate the model parameters and the acceleration factor. This method is more consistent and efficient, providing estimates with greater statistical precision and expressing uncertainty using confidence limits.
Assume that t 1,f , t 2,f , . . . , t r 1 ,f are the failure timings of r 1 units that have failed under normal condition X 1 , as well as m 1 censored units at periods t 1,s , t 2,s , . . . , t m 1 ,s . We further suppose that t 1,f , t 2,f , . . . , t r 2 ,f are the r 2 failure times at accelerated condition X 2 with m 2 censored units with censoring times t 1,s , t 2,s , . . . , t m 2 ,s . Now, for multiply censored data, the likelihood function under SSALT may be stated as follows Wang et al. 21 : The log-likelihood function l = L(t, α, β, θ ) based on MCS under SSPALT corresponding to Eq. (8) after substituting the values of f t i,f , f t k,f , F t j,s and F t l,s can be expressed as: where (τ + θ(−τ + t k )) = A k and (τ + θ(−τ + t l )) = A l . Now, the likelihood equations may be derived by calculating partial derivatives of Eq. (9) with respect to α, β and θ as: The ML estimates of the parameters, say α,β and θ , may be derived by solving Eq. (10)- (12) with regard to α, β and θ respectively. These equations are extremely complex, and they cannot be solved analytically. To solve these simultaneous equations, a numerical iteration approach such as a generic method named as Nelder-Mead Method is advised. In this paper, we utilized the Optim() function of R Software.  www.nature.com/scientificreports/ Because the precise sample distribution of the ML estimates cannot be determined in closed form, the estimated confidence intervals for the parameters α, β and θ are derived by utilizing the approximate distributions of their ML estimates, which is required to compute the Fisher information matrix. Because the predicted information matrix is excessively complex and necessitates numerical integration, the observed information matrix is produced. The asymptotic distribution of ML estimates of α, β and θ is given as α − α , β − β , θ − θ → N 0, I −1 (α, β, θ ) , where I represent 3 × 3 observed information matrix given in the following equation and the elements of I are given in the below.
And I −1 (α, β, θ ) represent the variance-covariance matrix and may be approximated by the inverse of the information matrix I, which contains the unknown parameters α, β and θ . Due to the consistency of the ML estimators, the parameters are substituted by the necessary MLEs to obtain an estimated I −1 (α, β, θ ) , which is given by Î −1 α,β,θ as follows: Now, the estimated 100(1 − ψ) percent double-sided confidence bounds with Z ψ/2 as an upper (ψ/2)th percentile of standard normal variate for α, β and θ are presented as follows: Second derivatives.

Real engineering application
In this section, we employ real data supplied by Proschon 36 to demonstrate the real engineering application of the estimation approaches offered in this paper. The data set provides the times of successive failures of sequential air conditioning system failures for each member of a Boeing 720 jet aircraft fleet. This data is also saved in R's npsurv package under the name acfail Wang 37 . The recorded data is given as follows: 1 1 2 3 3 3 3 4 5 5 5 5 5 7 7 7 9 9 10 11 11 11 11 12 12 12 12 13 14 14 14 14 14 14 14 14 15 15 15 16 16 16 18 18  18 18 18 18 20 20 21 21 22 22 22 23 23 23 24 24 25 26 26 27 27 29 29 29 29  To assess the data's goodness-of-fit to the NH distribution, we first calculated the NH distribution's parameters and then the K-S test was utilized. The K-S statistic and their p value are then calculated and reported in Table 1 as follows: The K-S distance is determined to be 0.04613, with an associated p value of 0.7552. Since the p value is more than 0.05, we cannot reject the null hypothesis, which states that both the theoretical and sample distributions www.nature.com/scientificreports/ are same. Furthermore, multiple plots are analyzed for the goodness-of-fit test to check further if the data fits the NH distribution. We also illustrate the fitting of the distribution by plotting the estimated cdf, Q-Q and P-P plots of the NH distribution for the supplied real data set. Figure 2 compares the theoretical CDF of the NH distribution to the empirical CDF and histogram. The Q-Q and P-P plots of the supplied actual data set are shown in Fig. 3. The figures and KS test results suggest that the real data set under consideration fits the NH distribution pretty well. Now, under SSPALT, we assume that the test runs under normal operating condition until the high stress is applied, and then the stress is raised to make the test operate at accelerated condition. We assume that the stress change time τ is 57 and various censoring levels (CL) are 20%, 30%, and 40%. To reflect multiple censoring strategy, we now delete 20%, 30%, and 40% of data at both stress levels, and the observed data under standard and accelerated stress with the assumed censoring levels are provided in Table 2.
We now estimated the MLEs and loglikelihood function values with stress change time τ = 57 and various censoring levels using the real data set reported in Table 2. The findings are computed using the R software, and Table 3 summarizes the resulting MLEs and loglikelihood function values. As per the results in Table 3, the estimates perform better at a 20 percent censoring level than at 30 percent and 40 percent filtering levels. This is apparent since larger data sets yield better results with more precision. We have now raised the stress change duration from 57 to 71 while keeping the same censoring levels of 20%, 30%, and 40% to test the model's flexibility. Table 4 shows the observed data under standard and accelerated stress, as well as the censoring levels after eliminating 20%, 30%, and 40% of the data for both stress levels.
Subsequently, we estimated the MLEs and loglikelihood function values with stress change time τ = 71 and various censoring levels using the real data set reported in Table 4. The findings are again computed using the R software, Team 38 and Table 5 summarizes the resulting MLEs and loglikelihood function values. As per the results in Table 5, we again observed that the estimates perform better at a 20% censoring level than at 30%, and 40% filtering levels.

Simulation study
This section provides a simulation study to investigate the performance of MLEs of parameters for NH distribution under SSPALT based on multiply censored data. The mean squared error (MSE) is used to compare the performance of point estimates, whereas the CPs are used to compare the performance of ACIs. The data were derived from Eq. (2) using the inverse CDF approach. The quantile function that has been utilized for this task is defined as t = 1 − log(1 − u) 1/α − 1 /β , where u is generated from uniform distribution, i.e., u ∼ U(0, 1) .
We selected five different sample sizes n = 80, 90, 100, 110, and 120 in order to analyse the nature of the estimates as the sample size increased. We also considered the NH distribution in order to produce data with starting values of the shape parameter α = 0.2 and the scale parameter β = 1.6 . The value of acceleration factor is set to θ = 2.5 , with two distinct values of stress change time τ = 5, 8 and three different censoring levels (CL) of 20%, 30%, and 40%. The simulation study is developed based on 8,000 multiply censored samples under SSPALT to observe changes in parameter values. The complete steps of the algorithm are detailed below: i. Using the pre-specified parameter values, generate n random samples from the NH distribution. To do so, first generate u from a uniform distribution with the command u = runif (n, 0, 1) , and then use the quantile function t = 1 − log(1 − u ) 1/α − 1 /β to generate the values of t = (t 1 , t 2 , . . . , t n ) from an NH distribution. ii. Now select the number of failures before the stress change time τ and denote it as n 1 . The generated sample up to time τ is t 1s = t 1 , t 2 , . . . , t n 1 . Also select the number of failed items at accelerated condition and denote it as n 2 , where n 1 + n 2 = n and the generated sample is given as t 2s = t n 1 +1 , t n 1 +2 . . . , t n 1 +n 2 Table 2. Multiply censored data at both stress level with τ = 57 and CLs of 20%, 30%, and 40%.

Conclusions and suggestions for further studies
In this article, we used TRV modeling for SSPALT to estimate the unknown model parameters of the NH distribution using the MLE technique. It has been discovered that the MLEs for all unknown parameters cannot be derived explicitly. As a result, we utilized R software to compute MLEs numerically using the Optim() function.
Using the observed Fisher information matrix, ACIs for the unknown parameters were also computed. An actual data set based on the timings of subsequent failures of sequential air conditioning system failures for each member of a Boeing 720 jet aircraft fleet was analyzed to illustrate the applicability of the different techniques. Finally, extensive simulation tests with various censoring mechanisms were carried out to evaluate the performance of Table 8. MLEs, MSEs, lower 95% ACI limit (L95%CL), upper 95% ACI limit (U95%CL) and 95% ACI coverage probability (95%ACICP) under multiply censored data with (α = 0.2, β = 1.6, θ = 2.5, CL = 0.40, τ = 5).  Table 9. MLEs, MSEs, lower 95% ACI limit (L95%CL), upper 95% ACI limit (U95%CL) and 95% ACI coverage probability (95%ACICP) under multiply censored data with (α = 0.2, β = 1.6, θ = 2.5, CL = 0.20, τ = 8). www.nature.com/scientificreports/ the estimate procedure. In particular, MSEs, ACIs, and corresponding average interval lengths were used as benchmarks. According to our numerical findings, the values of MSEs and average lengths drop as sample size increases. Furthermore, when the censoring level is reduced, the considered estimates of α, β and θ approach to their real values. As a future research, researchers may use rank set sampling to examine the NH distribution for hybrid censored data under SSPALT. A Bayesian analysis may be performed and compared with present study for the multiple censoring technique. Same has been added in the "Conclusion" section.  Table 11. MLEs, MSEs, lower 95% ACI limit (L95%CL), upper 95% ACI limit (U95%CL) and 95% ACI coverage probability (95%ACICP) under multiply censored data with (α = 0.2, β = 1.6, θ = 2.5, CL = 0.40, τ = 8).

Data availability
All data available in the paper with references.

Author contributions
All authors reviewed the manuscript and contributed equally to this paper.